#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: xjx
@time: 2023/10/28 14:40 
@file: logistics_mode_model.py
@project: nanchangproject
@describe: 物流模式识别模型
"""
import math
import numpy as np
import sys
import os
import django
from django_pandas.io import read_frame
from utils.DrivingMode.driving_mode_model import StatisticalResult
import pandas as pd
import transbigdata as tbd

sys.path.append('../../')
os.environ.setdefault('DJANGO_SETTINGS_MODULE', 'nanchangproject.settings')
django.setup()  # 很关键的配置，必须先setup在导入models类

from drivinginfo.models import LogisticsMode, BasicInfo, Drivinginfo
from django.core.exceptions import ObjectDoesNotExist
from django.db.models import Sum, Avg, Max, Min


class LogisticsModeBase(object):
    def __init__(self, radius_of_gyration, daily_travel_location, motion_entropy,
                 daily_driving_distance, daily_driving_time, logistics_mode: str = ''):
        self.radius_of_gyration = radius_of_gyration  # 回转半径
        self.daily_travel_location = daily_travel_location  # 日均出行地点
        self.motion_entropy = motion_entropy  # 活动熵
        self.daily_driving_distance = daily_driving_distance  # 日均行驶里程
        self.daily_driving_time = daily_driving_time  # 日均行驶时间
        self.logistics_mode = logistics_mode  # 物流模式

    def to_list(self):
        return [self.radius_of_gyration, self.daily_travel_location, self.motion_entropy,
                self.daily_driving_distance, self.daily_driving_time]


class LogisticsModeInstance(LogisticsModeBase):
    def __init__(self, radius_of_gyration, daily_travel_location, motion_entropy,
                 daily_driving_distance, daily_driving_time, car_number, vin):
        super(LogisticsModeInstance, self).__init__(radius_of_gyration, daily_travel_location, motion_entropy,
                                                    daily_driving_distance, daily_driving_time)
        self.car_number = car_number
        self.vin = vin


class LogisticsModeUtil(object):
    def __init__(self):
        self.LogisticsModeCluster1 = LogisticsModeBase(radius_of_gyration=1.24, daily_travel_location=7.39,
                                                       motion_entropy=1.26, daily_driving_time=25685.04,
                                                       daily_driving_distance=248.83,
                                                       logistics_mode='长距离往返型')
        self.LogisticsModeCluster2 = LogisticsModeBase(radius_of_gyration=0.41, daily_travel_location=3.57,
                                                       motion_entropy=0.51, daily_driving_time=14268.86,
                                                       daily_driving_distance=61.19,
                                                       logistics_mode='中距离定点型')
        self.LogisticsModeCluster3 = LogisticsModeBase(radius_of_gyration=0.23, daily_travel_location=8.82,
                                                       motion_entropy=1.14, daily_driving_time=25690.94,
                                                       daily_driving_distance=100.59,
                                                       logistics_mode='短距离多点型')

    # 提取od数据，计算回转半径、日均出行地点等数据，根据这些维度进行聚类分析
    def logistics_mode_recognition(self, instance: LogisticsModeInstance):
        instance_list = instance.to_list()
        cluster_list1 = self.LogisticsModeCluster1.to_list()
        cluster_list2 = self.LogisticsModeCluster2.to_list()
        cluster_list3 = self.LogisticsModeCluster3.to_list()
        cov = np.var(np.vstack([[instance_list], [cluster_list1], [cluster_list2], [cluster_list3]]), axis=0, ddof=1)
        print(f'方差为：{cov}')
        sum1 = 0
        sum2 = 0
        sum3 = 0
        for i in range(len(instance_list)):
            if abs(cov[i]) < 1e-6:
                continue
            sum1 += (instance_list[i] - cluster_list1[i]) ** 2 / cov[i]
            sum2 += (instance_list[i] - cluster_list2[i]) ** 2 / cov[i]
            sum3 += (instance_list[i] - cluster_list3[i]) ** 2 / cov[i]
        distance1 = math.sqrt(sum1)
        distance2 = math.sqrt(sum2)
        distance3 = math.sqrt(sum3)
        dis_list = [(distance1, self.LogisticsModeCluster1.logistics_mode),
                    (distance2, self.LogisticsModeCluster2.logistics_mode),
                    (distance3, self.LogisticsModeCluster3.logistics_mode)]
        dis_list = sorted(dis_list, key=lambda dis: dis[0])
        print(dis_list)
        instance.logistics_mode = dis_list[0][1]

    @staticmethod
    def logistics_mode_load_to_database(instance: LogisticsModeInstance):
        # 查询basic_info数据库中是否有该vin码的数据，若有则返回该数据对象，若无则创建
        try:
            basic_instance = BasicInfo.objects.get(vin=instance.vin)
        except ObjectDoesNotExist:
            basic_instance = BasicInfo.objects.create(
                vin=instance.vin,
                car_number=instance.car_number
            )
        # 查询logistics_mode数据库中是否有该basic_info的数据，若有则更新，若无则创建
        LogisticsMode.objects.update_or_create(
            basic_info=basic_instance,
            defaults={
                'logistics_mode': instance.logistics_mode,
                'radius_of_gyration': instance.radius_of_gyration,
                'daily_travel_location': instance.daily_travel_location,
                'motion_entropy': instance.motion_entropy,
                'daily_driving_distance': instance.daily_driving_distance,
                'daily_driving_time': instance.daily_driving_time,
            }
        )
        print(f'{instance.vin}车辆物流模式识别数据写入成功')


def calculate_logistics_mode(vin: str):
    """
    用于计算超速时间占比等信息
    @return: None
    """
    qs = Drivinginfo.objects.filter(vin=vin).order_by('data_time')  # 升序排列
    max_time = qs.aggregate(Max('data_time'))
    min_time = qs.aggregate(Min('data_time'))
    max_year = list(max_time.values())[0].year  # int类型，最近一年
    min_year = list(min_time.values())[0].year  # 最远一年
    year_iterator = min_year
    dataframe = read_frame(qs=qs)  # 将数据库对象转换为dataframe对象
    calculate_result = []
    while True:
        # 统计物流模式应当以年为单位还是以月为单位？
        for month in range(1, 13):
            # 筛选当前月份的数据,判断是否有数据，无则跳过当前循环
            monthly_data = dataframe[
                (dataframe['data_time'].dt.month == month) & (dataframe['data_time'].dt.year == year_iterator)]
            if monthly_data.empty:
                continue
            df = monthly_data[
                ['vin', 'data_time', 'longitude', 'latitude', 'air_intake', 'engine_status',
                 'accumulation_mile']].copy()
            # 用进气量进行发动机状态缺失值补偿
            df.loc[(df['air_intake'] == '0.0') & (df['engine_status'] == 'nan'), 'engine_status'] = '关闭'
            df.loc[(df['air_intake'] != '0.0') & (df['engine_status'] == 'nan'), 'engine_status'] = '启动'
            level_map = {"启动": 1, "关闭": 0}
            df['engine_status'] = df['engine_status'].map(level_map)
            # 对定位缺失值进行补偿
            df['longitude'] = df['longitude'].replace(0, np.nan)
            df['latitude'] = df['latitude'].replace(0, np.nan)
            df['longitude'].fillna(method='bfill', inplace=True)
            df['latitude'].fillna(method='bfill', inplace=True)
            df['longitude'] = df['longitude'].astype(float)
            df['latitude'] = df['latitude'].astype(float)
            # 对累计里程缺失值进行补偿
            df['accumulation_mile'].fillna(method='bfill', inplace=True)
            df['data_time'] = pd.PeriodIndex(df['data_time'], freq="S")  # 以秒为间隔创建序列，可以自由的进行加减
            time_sum = (df['data_time'].iloc[0] - df['data_time'].iloc[-1]).n  # 该月时间间隔的总秒数
            # 计算回转半径
            # 计算 "经度 "和 "纬度"（质心坐标）的平均值
            x_bar = df['longitude'].mean()
            y_bar = df['latitude'].mean()
            # 计算每个数据点到质量中心的距离的平方
            distance_squared = (df['longitude'] - x_bar) ** 2 + (df['latitude'] - y_bar) ** 2
            # 计算回转半径（惯性矩）：距离平方的平均值并取平方根
            radius_of_gyration = round(np.sqrt(distance_squared.mean()), 6)  # 保留六位小数

            # 计算日均活动位置数
            # 提取OD点数据
            engine_status_df = df['engine_status'].copy()
            start_p = engine_status_df.diff()  # diff是当前数据减去上一个数据，第一个则是NaN
            index_start = df[start_p == 1].copy()
            index_stop = df[start_p == -1].copy()
            print(f'{month}月共有停车点{len(index_stop)}')
            # index_start和index_stop有可能是空，即当天没有停车点，此时应跳出循环
            if index_start.empty or index_stop.empty:
                continue

            # 计算日均行驶里程
            daily_driving_distance = round((df['accumulation_mile'].iloc[-1] - df['accumulation_mile'].iloc[0]) / 30, 6)

            # 筛选OD点数据
            index_start = index_start.reset_index(drop=True)
            index_stop = index_stop.reset_index(drop=True)
            if index_start['data_time'].iloc[0] > index_stop['data_time'].iloc[0]:
                # 注意顺序，应当从第一个启动点开始计算，因此第一个终止点应当大于第一个起始点
                index_stop = index_stop.drop(index=0, inplace=False)
                # 若不满足上述条件，应当删除第一个终止点,inplace=False表示返回一个copy
                index_stop = index_stop.reset_index(drop=True)  # 重置索引
            time_stay = index_stop['data_time'] - index_start['data_time']
            index_stop['duration_time'] = time_stay.copy()  # 行驶持续时间
            index_stop = index_stop.dropna(axis=0)  # 删除NaN的行
            index_stop['duration_time'] = index_stop['duration_time'].apply(lambda x: x.n)  # .n表示取秒数

            # 计算日均行驶时间
            daily_driving_time = index_stop['duration_time'].sum()  # 对行驶时间（秒）求和
            daily_driving_time = round(daily_driving_time / 30, 6)  # 除以30表示日均

            # 计算日均活动数
            selected_stop = index_stop[index_stop['duration_time'] >= 600]  # 行驶时间大于10分钟作为活动地点移动的阈值
            daily_travel_location = round(len(selected_stop) / 30, 6)  # 日均活动位置数

            # 使用transbigdata库进行地图栅格化
            data1 = df[(df['longitude'] != 0) & (df['latitude'] != 0)].copy()  # 去除异常定位数据
            data1['engine_status'] = data1['engine_status'].replace({0: 1, 1: 0})
            data1.columns = ['VehicleNum', 'Time', 'Lng', 'Lat', 'air_intake', 'OpenStatus', 'accumulation_mile']
            oddata = tbd.taxigps_to_od(data1, col=['VehicleNum', 'Time', 'Lng', 'Lat', 'OpenStatus'])
            oddata['duration_time'] = oddata['etime'] - oddata['stime']
            oddata['duration_time'] = oddata['duration_time'].apply(lambda x: x.n)
            # 计算活动边界
            k = oddata[oddata['duration_time'] >= 600]  # 行驶时间大于10分钟作为活动地点移动的阈值
            lng_min = df['longitude'].min()
            lat_min = df['latitude'].min()
            lng_max = df['longitude'].max()
            lat_max = df['latitude'].max()
            # 定义研究范围边界
            bounds = [lng_min, lat_min, lng_max, lat_max]
            # 通过边界获取栅格化参数
            params = tbd.area_to_params(bounds, accuracy=300)
            # 将GPS数据对应至栅格，将生成的栅格编号列赋值到数据表上作为新的两列
            data = pd.DataFrame()
            data['LONCOL'], data['LATCOL'] = tbd.GPS_to_grid(k['elon'], k['elat'], params)
            k.reset_index(drop=True, inplace=True)
            k = pd.concat([k, data], axis=1)
            # 计算重复位置点
            rep_point = data.duplicated()
            for i in range(len(k)):
                if rep_point.iloc[i] == True:
                    tmp_lng, tmp_lat = data['LONCOL'].iloc[i], data['LATCOL'].iloc[i]
                    reptmp = (k['LONCOL'] == tmp_lng) & (k['LATCOL'] == tmp_lat)
                    c = k.loc[reptmp, 'duration_time']
                    c.iloc[0] += k['duration_time'].iloc[i]
                    tmp = c.iloc[0]
                    k.loc[c.index[0], 'duration_time'] = tmp
            # 计算pi
            m_pi = k[rep_point == False].copy()
            pi_sum = m_pi['duration_time'].sum()
            dur_pi = m_pi['duration_time'].values
            dur_pi = dur_pi / pi_sum
            # 计算motion_entropy即活动熵
            motion_entropy = round(-1 * np.sum(dur_pi * np.log(dur_pi)), 6)

            if not np.isnan(daily_travel_location):
                calculate_result.append(StatisticalResult(year=year_iterator,
                                                          month=month,
                                                          what_data='daily_travel_location',
                                                          data=daily_travel_location))
            if not np.isnan(daily_driving_distance):
                calculate_result.append(StatisticalResult(year=year_iterator,
                                                          month=month,
                                                          what_data='daily_driving_distance',
                                                          data=daily_driving_distance))
            if not np.isnan(daily_driving_time):
                calculate_result.append(StatisticalResult(year=year_iterator,
                                                          month=month,
                                                          what_data='daily_driving_time',
                                                          data=daily_driving_time))
            if not np.isnan(motion_entropy):
                calculate_result.append(StatisticalResult(year=year_iterator,
                                                          month=month,
                                                          what_data='motion_entropy',
                                                          data=motion_entropy))
            if not np.isnan(radius_of_gyration):
                calculate_result.append(StatisticalResult(year=year_iterator,
                                                          month=month,
                                                          what_data='radius_of_gyration',
                                                          data=radius_of_gyration))

        year_iterator += 1
        if year_iterator > max_year:
            break
    return calculate_result


def calculate_logistics_mode_all(vin: str):
    """
    以所有数据计算车辆的物流模式
    @param vin: 该车辆VIN码
    @return: 一个字典
    """
    qs = Drivinginfo.objects.filter(vin=vin).order_by('data_time')  # 升序排列
    dataframe = read_frame(qs=qs)  # 将数据库对象转换为dataframe对象
    if dataframe.empty:
        return
    df = dataframe[
        ['vin', 'data_time', 'longitude', 'latitude', 'air_intake', 'engine_status',
         'accumulation_mile']].copy()
    # 用进气量进行发动机状态缺失值补偿
    df.loc[(df['air_intake'] == '0.0') & (df['engine_status'] == 'nan'), 'engine_status'] = '关闭'
    df.loc[(df['air_intake'] != '0.0') & (df['engine_status'] == 'nan'), 'engine_status'] = '启动'
    level_map = {"启动": 1, "关闭": 0}
    df['engine_status'] = df['engine_status'].map(level_map)
    # 对定位缺失值进行补偿
    df['longitude'] = df['longitude'].replace(0, np.nan)
    df['latitude'] = df['latitude'].replace(0, np.nan)
    df['longitude'].fillna(method='bfill', inplace=True)
    df['latitude'].fillna(method='bfill', inplace=True)
    df['longitude'] = df['longitude'].astype(float)
    df['latitude'] = df['latitude'].astype(float)
    # 对累计里程缺失值进行补偿
    df['accumulation_mile'].fillna(method='bfill', inplace=True)
    df['data_time'] = pd.PeriodIndex(df['data_time'], freq="S")  # 以秒为间隔创建序列，可以自由的进行加减
    time_sum = (df['data_time'].iloc[0] - df['data_time'].iloc[-1]).n  # 该月时间间隔的总秒数
    # 计算回转半径
    # 计算 "经度 "和 "纬度"（质心坐标）的平均值
    x_bar = df['longitude'].mean()
    y_bar = df['latitude'].mean()
    # 计算每个数据点到质量中心的距离的平方
    distance_squared = (df['longitude'] - x_bar) ** 2 + (df['latitude'] - y_bar) ** 2
    # 计算回转半径（惯性矩）：距离平方的平均值并取平方根
    radius_of_gyration = round(np.sqrt(distance_squared.mean()), 6)  # 保留六位小数

    # 计算日均活动位置数
    # 提取OD点数据
    engine_status_df = df['engine_status'].copy()
    start_p = engine_status_df.diff()  # diff是当前数据减去上一个数据，第一个则是NaN
    index_start = df[start_p == 1].copy()
    index_stop = df[start_p == -1].copy()
    print(f'共有停车点{len(index_stop)}')
    # index_start和index_stop有可能是空，即当天没有停车点，此时应跳出循环
    if index_start.empty or index_stop.empty:
        return
        # 计算日均行驶里程
    daily_driving_distance = round((df['accumulation_mile'].iloc[-1] - df['accumulation_mile'].iloc[0]) / 30, 6)

    # 筛选OD点数据
    index_start = index_start.reset_index(drop=True)
    index_stop = index_stop.reset_index(drop=True)
    if index_start['data_time'].iloc[0] > index_stop['data_time'].iloc[0]:
        # 注意顺序，应当从第一个启动点开始计算，因此第一个终止点应当大于第一个起始点
        index_stop = index_stop.drop(index=0, inplace=False)
        # 若不满足上述条件，应当删除第一个终止点,inplace=False表示返回一个copy
        index_stop = index_stop.reset_index(drop=True)  # 重置索引
    time_stay = index_stop['data_time'] - index_start['data_time']
    index_stop['duration_time'] = time_stay.copy()  # 行驶持续时间
    index_stop = index_stop.dropna(axis=0)  # 删除NaN的行
    index_stop['duration_time'] = index_stop['duration_time'].apply(lambda x: x.n)  # .n表示取秒数

    # 计算日均行驶时间
    daily_driving_time = index_stop['duration_time'].sum()  # 对行驶时间（秒）求和
    daily_driving_time = round(daily_driving_time / 30, 6)  # 除以30表示日均

    # 计算日均活动数
    selected_stop = index_stop[index_stop['duration_time'] >= 600]  # 行驶时间大于10分钟作为活动地点移动的阈值
    daily_travel_location = round(len(selected_stop) / 30, 6)  # 日均活动位置数

    # 使用transbigdata库进行地图栅格化
    data1 = df[(df['longitude'] != 0) & (df['latitude'] != 0)].copy()  # 去除异常定位数据
    data1['engine_status'] = data1['engine_status'].replace({0: 1, 1: 0})
    data1.columns = ['VehicleNum', 'Time', 'Lng', 'Lat', 'air_intake', 'OpenStatus', 'accumulation_mile']
    oddata = tbd.taxigps_to_od(data1, col=['VehicleNum', 'Time', 'Lng', 'Lat', 'OpenStatus'])
    oddata['duration_time'] = oddata['etime'] - oddata['stime']
    oddata['duration_time'] = oddata['duration_time'].apply(lambda x: x.n)
    # 计算活动边界
    k = oddata[oddata['duration_time'] >= 600]  # 行驶时间大于10分钟作为活动地点移动的阈值
    lng_min = df['longitude'].min()
    lat_min = df['latitude'].min()
    lng_max = df['longitude'].max()
    lat_max = df['latitude'].max()
    # 定义研究范围边界
    bounds = [lng_min, lat_min, lng_max, lat_max]
    # 通过边界获取栅格化参数
    params = tbd.area_to_params(bounds, accuracy=300)
    # 将GPS数据对应至栅格，将生成的栅格编号列赋值到数据表上作为新的两列
    data = pd.DataFrame()
    data['LONCOL'], data['LATCOL'] = tbd.GPS_to_grid(k['elon'], k['elat'], params)
    k.reset_index(drop=True, inplace=True)
    k = pd.concat([k, data], axis=1)
    # 计算重复位置点
    rep_point = data.duplicated()
    for i in range(len(k)):
        if rep_point.iloc[i] == True:
            tmp_lng, tmp_lat = data['LONCOL'].iloc[i], data['LATCOL'].iloc[i]
            reptmp = (k['LONCOL'] == tmp_lng) & (k['LATCOL'] == tmp_lat)
            c = k.loc[reptmp, 'duration_time']
            c.iloc[0] += k['duration_time'].iloc[i]
            tmp = c.iloc[0]
            k.loc[c.index[0], 'duration_time'] = tmp
    # 计算pi
    m_pi = k[rep_point == False].copy()
    pi_sum = m_pi['duration_time'].sum()
    dur_pi = m_pi['duration_time'].values
    dur_pi = dur_pi / pi_sum
    # 计算motion_entropy即活动熵
    motion_entropy = round(-1 * np.sum(dur_pi * np.log(dur_pi)), 6)
    return {
        'radius_of_gyration': radius_of_gyration,
        'daily_travel_location': daily_travel_location,
        'motion_entropy': motion_entropy,
        'daily_driving_distance': daily_driving_distance,
        'daily_driving_time': daily_driving_time,
    }

def generate_logistics_mode(vin: str):
    """
    用于生成或更新车辆的物流模式识别结果
    @param vin:
    @return:
    """
    result = calculate_logistics_mode_all(vin=vin)
    print(result)
    instance = LogisticsModeInstance(
        radius_of_gyration=result['radius_of_gyration'],
        daily_travel_location=result['daily_travel_location'],
        motion_entropy=result['motion_entropy'],
        daily_driving_distance=result['daily_driving_distance'],
        daily_driving_time=result['daily_driving_time'],
        car_number='',
        vin=vin
    )
    x = LogisticsModeUtil()
    x.logistics_mode_recognition(instance=instance)
    print(instance.logistics_mode)
    x.logistics_mode_load_to_database(instance=instance)
    print(f'车辆{vin}驾驶风格识别成功，结果{instance.logistics_mode}')


if __name__ == '__main__':
    generate_logistics_mode(vin='LETYBEG27LH039366')



